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Abstract. We calculate the two-, three- and (for the first time) four-point correlation functions of the COBE-DMR 
4-year sky maps, and search for evidence of non-Gaussianity by comparing the data to Monte Carlo-simulations of 
the functions. The analysis is performed for the 53 and 90 GHz channels, and five linear combinations thereof. For 
each map, we simulate an ensemble of 10 000 Gaussian realizations based on an a priori best-fit scale-invariant 
cosmological power spectrum, the DMR beam pattern and instrument-specific noise properties. Each observed 
COBE-DMR map is compared to the ensemble using a simple \ 2 statistic, itself calibrated by simulations. In 
addition, under the assumption of Gaussian fluctuations, we find explicit expressions for the expected values of the 
four-point functions in terms of combinations of products of the two-point functions, then compare the observed 
four-point statistics to those predicted by the observed two-point function, using a redefined \ 2 statistic. Both 
tests accept the hypothesis that the DMR maps are consistent with Gaussian initial perturbations. 
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1. Introduction 

The study of CMB temperature anisotropies and their sta- 
tistical properties has become an important theme in mod- 
ern cosmology. In its most conventional interpretation, the 
distribution of anisotropies reflects the properties of the 
universe approximately 300 000 years after the Big Bang, 
at the surface of last scattering. Thus, by measuring sta- 
tistical quantities such as the angular power spectrum or 
the angular two-point correlation function, we can infer 
the values of many interesting cosmological parameters. 

For both theoretical and practical purposes, it is con- 
venient to expand the temperature anisotropy field into a 
sum of (complex) spherical harmonics: 
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The temperature perturbation field is said to be Gaussian 
distributed if each ai m follows an independent Gaussian 
probability distribution. The question of whether the ob- 
served temperature field is Gaussian or otherwise is of 
crucial importance for modern cosmology. From a scien- 
tific standpoint, most conventional inflationary models of 
structure formation predict a Gaussian temperature field, 
whereas scenarios which invoke topological defects to seed 
the large-scale structure predict a non-Gaussian distribu- 
tion. Thus the statistical properties of the a; m 's can be 



used to distinguish such models. Secondly, from an anal- 
ysis point-of-view, most parameter estimation techniques 

- usually based on the observed angular power spectrum 

- assume Gaussianity, and may therefore be biased if the 
observed field is indeed non-Gaussian. 

However, testing for non-Gaussianity is anything but 
trivial, and several qualitatively different tests are re- 
quired in order to perform a complete analysis. At present 
the arsenal of available tests which have been applied to 
the COBE-DMR data consists of at least the follo wing : 
bi- and tr ispec trum based analysis (F erreir a et al. 1998 ; 
Ma gueijo 2000 : Sand vik fc Magueijo 2001 ; Komatsu et 
al. 2002; Kunz et al. 2001), 3-point correlation function 



based tests (Kogut et al. 1996), methods utilizing wavelets 
(Cayon et al. 2001; Barreiro et al. 
functionals (Schmalzing & Gorski 



200C) 



and Minkowski 
19981; Novikov et al. 



[2000 ) . Indeed, there has been a small resurgence in interest 
in the possibility of non-Gaussian signals in the COBE- 
DMR maps as a consequence of the bispectrum work of 
Ferreira et al. (1998) and Magueijo ( pOOOl ). These papers 
find non-Gaussian contributions using harmonic analyses 
at the 98% confidence limit, and although Banday et al. 
( 2000D explain these tentative detections by appealing to 
the presence of a specific residual systematic artifact in 
the data, additional investigation is warranted. 
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In this paper we adopt A-point correlation functions 
as probes for non-Gaussianity. For a Gaussian field all 
odd TV-point functions (such as the three-point function) 
have vanishing expectation values, while all even A-point 
functions can be reduced to expressions involving the two- 
point function. Thus, if the observed three-point function 
is significantly non-zero when compared to a Gaussian en- 
semble, its native distribution is probably non-Gaussian. 
Further, if the four-point function does not reduce into 
two-point functions, the same conclusion can be made. 

The first part of this paper builds on ideas demon- 
strated in Kogut et al. ( |1996| ) and Hinshaw et al. ( 1995| ). 
We study the 4-year COBE-DMR sky maps, computing 
the two- and three-point functions, as has been performed 
previously, then proceeding to extend the analysis for the 
first time to the determination of several four-point func- 
tions. The definitions of these new functions are given in 
Sect. H 

In Sect. |^ we compute the various correlation functions 
for the four DMR channels and five linear combinations 
thereof. Next, we compute the same functions for 10 000 
Monte Carlo simulated Gaussian maps, which are used as 
the basis of the various statistical tests of non-Gaussianity. 
Initially, we apply the x 2 test as defined by Kogut et al. 



1996, by comparing the observed data value to the distri- 



bution generated from the application of the x statistic 
for each map in the simulated ensemble. 

Subsequently, in Sect. |j, we provide expressions for the 
expected value of several four-point functions in terms of 
the two-point function, then explicitly compare the ob- 
served four-point functions to those predicted by the ob- 
served two-point function. We define a suitable % 2 statistic 
in order to quantitatively measure the degree of deviation, 
once again calibrated by Monte Carlo-simulations. 

2. Definitions of the correlation functions 

An A-point correlation function is defined as the average 
product of A temperatures with a fixed relative orienta- 
tion on the sky: 
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where hi is the direction vector of the z'th pixel in the 
configuration, and cos 8j = hk • ni for 27V — 3 arbitrary, 
but different, angular pixel distances on the sky. 

Although the TV-point functions are easily defined 
and relatively simple to implement computationally, their 
evaluation is generally CPU-intensive, which is especially 
problematic since detailed assessment of results requires 
large Monte Carlo simulation data sets. The full computa- 
tion of an TV-point correlation function scales as 0(N^ 1X ), 
and is therefore virtually impossible to compute for high- 
resolution maps for any order A greater than two. For this 
reason we choose to compute only a subset of the possibili- 
ties in the TV-dimensional configuration space, designed to 
reduce the complexity of the problem. As an example con- 
sider the pseudo-collapsed three-point function for which 



we require two points to coincide, effectively reducing the 
geometry to that of the two-point function. Such subsets 
typically scale somewhere between 0(Ap ix ) and 0(Ap ix ). 
Thus, with some effort put into the implementation these 
functions can be computed even for rather high-resolution 
maps. 

Previous work has considered two special three-point 
functions, namely the collapsed and the equilateral func- 



tions (Kogut et al. 1996; Hinshaw et al. 1995). As men- 
tioned above, the collapsed function is defined by requir- 
ing two of the three point to coincide, while the equilateral 
function requires the three points to span an equilateral 
triangle on the sphere. 

In this paper, we shall also consider several simple four- 
point configurations. These functions are, in order of com- 
plexity: 

1. the collapsed 1+3 point function 

2. the collapsed 2+2 point function 

3. the collapsed equilateral four-point function 

4. the rhombic four-point function 

The names should be self-explanatory. For the 1+3 point 
function, three points coincide in a manner similar to 
the collapsed three-point function. The 2+2 point func- 
tion is defined by allowing two pixel pairs to coincide. 
The collapsed equilateral function is the equilateral three- 
point configuration where one pixel is multiplied twice. 
The last case, the rhombic four-point function, consists of 
two equilateral triangles "glued" together along one side, 
effectively spanning a rhombus on the sphere. Note that 
these functions are chosen because of ease of implementa- 
tion, not because they are better suited for the testing of 
Gaussianity than other configurations. 

Several of the functions defined above are so-called col- 
lapsed functions, i.e. one pixel is multiplied one or more 
times with itself. Unfortunately, for noisy maps this ren- 
ders the function completely noise dominated. To rem- 
edy this problem we substitute the collapsed functions 
by so-called pse udo-c ollapsed versions, as introduced by 
Hinshaw et al. (|1995|) . For the COBE-DMR experiment 
the beam size is approximately 7°, while the pixel size 
is - necessarily for adequate sampling - ~ 1.8° (for the 
HEALPix A S ido = 32 pixelization used here). Therefore 
the CMB signal component between two neighboring pix- 
els is highly coherent, whereas the noise contributions are 
independent. Thus, instead of multiplying a given pixel 
by itself several times, we multiply the pixel by one or 
more of its immediate neighbors, then sum over all such 
possible products, effectively multiplying by an average 
over the nearest neighbors. Hence, we more generally de- 
fine a pseudo-collapsed function as an average product of 
pixels where at least one pixel is multiplied in the pseudo- 
collapsed sense, ie. by an average over its neighbors. The 
golden rule for our analysis is that no pixel is ever mul- 
tiplied with itself. This definition is then not completely 
equivalent to that introduced by Hinshaw et al. ( 1995[ ) 
They defined the pseudo-collapsed function as the average 
product of 1) a center pixel, 2) one of its neighbors and 3) 
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a far point, where the far point was not allowed to be the 
center pixel. However, it was allowed to be the neighboring 
pixel. Although not a major problem for the three-point 
function, we have determined that the inclusion of such 
a product renders the first bin of the four-point functions 
completely noise dominated. 

We also introduce one further small change compared 
to Hinshaw et al. (1995) in that we exclude the zeroth 
angular bin (for which all A pixels coincide) as any cos- 
mological information here is heavily suppressed due to 
the low signal-to-noise ratio. Indeed, the inclusion of the 
zeroth bin only acts to increase the variance of the \ 2 
statistic, and is therefore better omitted. 



3. Measurements of the COBE-DMR correlation 
functions 

The COBE-DMR experiment resulted in six independent 
maps, two for each of the three frequencies at 31.5, 53 and 
90 GHz. In this work we only include the maps from the 
53 and 90 GHz channels, as they are superior in terms 
of the signal-to-noise ratio. The maps are analyzed in the 
HEALPixQ pixelization scheme, with a resolution param- 
eter of Aside = 32, corresponding to 12 288 pixels on the 
sky. At each frequency we compute the 'sum' (A + B)/2 
and 'difference' (A — B)/2 combinations, which yield, re- 
spectively, maps with enhanced signal-to-noise or noise 
content alone. In addition, we also generate a co-added 
map from the four basic channels using weights to achieve 
the optimal signal-to-noise ratio. 

The two-, three- and four-point correlation functions 
for these nine map combinations are then computed. All 
pixels corresponding to the extended Galactic cut (Banday 
et al. 



1997 but recomputed explicitly for the HEALPix 



scheme) are rejected from the analysis, leaving a total of 
7880 accepted pixels. The best-fitting monopole, dipole 
and quadrupole are subtracted from each map before the 
A-point functions are evaluated. The observed correlation 
functions are also computed after correction for the diffuse 
foreground emission at high Galactic latitude, using infor- 
mation from the ( appro priately scaled) DIRBE 140 /im 
map (Gorski et al. 1996 ). 

For our Monte Carlo ensemble, we simulate 10 000 in- 
dividual realizations of the CMB sky, based on an a pri- 
ori best-fit cosmological power spectrum. In particular, 
we consider scale-invariant Gaussian temperature fluctu- 



ations (P(fc) oc Q, 2 p R 
18 /aK, (Gorski et al. 1996). The power-spectrum is filtered 
through the DMR beam and pixel window functions. To 
each simulated CMB sky, we add four noise realizations 
based on the rms noise levels and observation patterns of 
the observed 53 and 90 GHz sky maps. These are then 
combined to generate the corresponding sum, difference 
and co-added sky maps. These are then processed in an 
identical fashion to the DMR data. 



k n with n = 1) with Q rms -ps 



We note that we have also assumed that there are no 
significant pixel-pixel noise correlations, although some 
will indeed be present as a consequence of the differen- 
tial nature of the radiometers which couple observations 
separated by ~ 60° on the sky. Lineweaver et al. ( 1994 ) 
have investigated this effect in detail, and find a small ex- 
cess signal is present in the 2-point correlation function 
at 60° for maps containing noise signal alone. However, 
we do not expect our results to be compromised by this 
assumption. 

Fig. [j] shows the results from these calculations for the 
co-added maps. The observed functions lie comfortably 
within the confidence region defined by the Monte Carlo- 
simulations and there are no striking deviations visible by 
simple inspection. 



4. The x 2 statistic 

In order to quantitatively measure the agreement between 
the DMR maps and the simulated ensemble, we utilize the 



same x 2 methodology described by Kogut et al. (1996): 



X 2 - £(A» - (5 a ))(M- 1 ) Q0 ( J D /3 - (^» (3) 

a/3 

Here a and /3 denote angular bins, D the DMR correla- 
tion function, and (So) the mean function computed from 
simulations. M is the binned covariancc matrix: 



(4) 



http://www. eso. org/ 'science /healpix/ 



The probability density function of this \ 2 statistic is 
established by computing the same statistic for all maps in 
the ensemble, (ie. by substituting D with each of the sim- 
ulated maps, S l ). The resulting histogram represents the 
probability distribution function against which we com- 
pare the values from the DMR maps. Table [j] records the 
fraction of simulated maps with higher x 2 values than the 
observed DMR map. Thus, values of order 0.01 or 0.99 
can be considered suspicious, while anything from 0.05 to 
0.95 is acceptable. 

In Kogut et al. ( |1996| ) the results for the pseudo- 
collapsed and the equilateral three-point functions are 
given for the 53 GHz (A + B)/2 map; they find the frac- 
tions to be respectively 0.66 and 0.31, while we find 0.65 
and 0.29. Considering the minor changes in the definitions 
of the correlation functions and the different pixelizations 
used, the agreement is most satisfactory. 

Overall, the numbers indicate that the COBE-DMK 
maps agree very well with the simulations. The optimal 
co-added map, for which the signal-to-noise ratio is the 
highest, returns results comfortably in the accepted range, 
as does the combined analysis of all A-point functions. 
We conclude that the DMR maps are compatible with the 
Gaussian hypothesis as measured by this test. 
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Fig. 1. Three- and four-point correlation functions of the co-added 4-year DMR map. Solid line shows the most likely 
value for each bin, dark shading shows the 68% confidence region and light shading the 95% confidence region, as 
computed by Monte Carlo-simulations. Dots represent functions for the uncorrected co-added map, and boxes shows 
the functions for the map for which high-latitude Galactic emission has been removed. Note the different angular units 
on the horizontal axis, reflecting the fact that the various functions are defined on different angular intervals. 
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Table 1. Results from \ 2 tests. The numbers indicate the fraction of simulated realizations with \ 2 value higher than 
for the respective COBE map. The lower half shows the results for the DMR maps after correction for high latitude 
Galactic emission. The upper half shows the results for the uncorrected maps. The effect of Galactic emission appears 
to be minimal: this is not unexpected since a best-fit quadrupole has been removed from the data before analysis, and 
the Galactic emission is dominated by such large-scale structure. 



5. Reducing four-point functions into two-point 
functions 

Since it may be noted that all even-ordered iV-point func- 
tions have non-vanishing expectation values determined 
by the two-point function, we can define an additional test 
of Gaussianity for the four-point functions. Explicitly, we 
take advantage of the following property (as, for exam- 
ple, described in Adler 1981 ): If Y\, Y 2 , ■ ■ ■ ,Y n is a set of 
real-valued random variables having a joint Gaussian dis- 
tribution and zero means, then for any integer m: 



{YiY 

(YxY 2 ---Y 2m ) = 



E 



' ^2m+l) 

(Y h Y i2 ) 







(5) 
(6) 



Here the sum goes over all (2m — 1)!! different ways of 
grouping Y±, Y 2 , . . . , Y 2m into m pairs. 

Thus, if the CMB temperature anisotropy field is in 
fact Gaussian distributed with zero mean, then all even 
iV-point functions can be reduced to combinations of prod- 
ucts of the two-point function. In particular, the four-point 
function reduces to: 



(T1T2T3T4) = 



(T 1 r 2 )<T 3 r 4 ) 

(TiT 3 )(T 2 T 4 ) 

(r 1 r 4 )(T 1 2 r 3 ) 



(7) 



For the four functions we defined in Sect. |^, we find the 
following expressions: 

C 1+3 {9) = 3C 2 (0)C 2 (6) 
C 2+2 {6) = C 2 (Q) 2 + 2 C 2 (6) 2 
C CCiui (6) = C 2 (0)C 2 (6) + 2C 2 (6) 2 
C rhomb (60 = C 2 {6) C 2 {6') + 2 C 2 {6f 



(8) 
(9) 
(10) 

(11) 



In equation ( |TT| ) 9' denotes the length of the longest axis 
of the rhombus, which is easily computed by spherical 
trigonometry: 

6' cos9 

cos T = e ( 12 ) 

2 cos § 

Note that the functions given by Eqs. (^-([Ti"|) should 
be interpreted as expectation values. That is, the observed 
four-point function should equal that given by Eq. (0) to 
the same extent that the three-point function equals zero. 
Thus, to establish an acceptable distribution for these rela- 
tions we again utilize our Monte-Carlo simulation set. We 
compare the simulated ensemble of functions to those pre- 
dicted by Eqs. (p[)-(pd|): for each realization in the ensem- 
ble, we compute the four-point function predicted by the 
two-point function relation above, and then evaluate the 
difference between this predicted and the observed four- 
point function. From these differences we generate 68% 
and 95% confidence intervals. Finally, the procedure is re- 
peated for the DMR maps, and the results are compared 
to the derived confidence intervals. 

This procedure has one major advantage compared to 
the one described in Sect. |3| the power spectrum only 
mildly affects the result. That is, the two most important 
contributions to the analysis come from the map itself, in 
the form of a two-point and a four-point function. The 
assumed power spectrum is only used for estimating the 
acceptable deviations, not the overall shape. Therefore, 
this procedure provides a more direct test for Gaussianity 
than the previous one. 

The results are shown for the co-added map in Fig. ^. 
The observed function lies well within the confidence re- 
gions about the predicted function for all four cases. For a 
more quantitative measure of the perceived agreement, we 
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Fig. 2. Comparison of the observed and the "reduced" four-point functions for the co-added DMR sky map. Boxes 
indicate the function predicted by the two-point functions, while shaded areas represent the 68% and 95% confidence 
intervals computed from Monte Carlo-simulations. The observed four-point functions are shown with a solid line. 
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Table 2. Results for the modified \ 2 statistic. The numbers refer to the fraction of simulations with \ 2 values higher 
than the corresponding COBE-DMR maps. The upper half shows the results for the uncorrected maps, while Galactic 
emission has been corrected for in the lower half. 



define a \ 2 statistic, incorporating the new degree of free- new function: 

dom provided by the predicted four-point function by sim- 2 _ V^m D prcd )fM _1 ) (D £)P rcd ) (13) 

ply replacing the average correlation function with that ' a a a P P p 
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where 



M, 



a/3 



^E(^-^ prcd )(^-^ prod ) 



(14) 



The meaning of each symbol is the same as in Eqs. (|[) 
and (|). 

Table |^ summarizes the results, which again support 
the hypothesis that the DMR sky maps are consistent with 
a scale-invariant cosmological model with Gaussian initial 
fluctuations. 

6. Conclusions 



By performing Monte Carlo-simulations we have studied 
the statistical properties of the COBE-DMR 53 GHz and 
90 GHz channels. The basic ingredients for this analysis 
were various iV-point correlation functions, and, in partic- 
ular, four different four-point functions which have been 
presented for the first time. We have additionally taken ad- 
vantage of a result from statistical theory, relating all even 
iV-point functions to reductions in terms of the two-point 
function. This allowed us to define a test for Gaussianity 
in which the assumed power spectrum only plays a sec- 
ondary role. This test could therefore prove better suited 
for situations in which we do not have access to the opti- 
mal power spectrum. 

Comparison of the DMR TV-point correlation func- 
tions with the Monte-Carlo ensemble indicates no evidence 
for possible non-Gaussian behavior, in agreement with 
the earlier analysis of Kogut et al. (1996). Furthermore, 
the agreement between the observed DMR functions and 
the simulated ensembles also supports the validity of our 
model assumptions, namely that of a scale-invariant power 
law model for the anisotropies, and uncorrelated noise. 

On the other hand, the excellent agreement between 
the simulated and the observed correlation functions poses 
an intriguing problem: tests of Gaussianity based on a har- 
monic analysis of th e DMR data - the bispectrum work of 
Fcrrei ra et al. ( |l998|) and trispectrum results of Kunz et al. 
( 2001 ) - show compelling evidence for non-Gaussian fea- 
tures (although these have subsequently been associated 
wi th sys tematic artifacts in the DMR data by Banday et 
al. 200C| ), while tests based on real-space high-order statis- 
tics such as those presented here do not. The resolution of 
such apparently contradictory results is most likely rather 
mundane: the source of the non-Gaussian signal was found 
to be strongly located at the multipole order I = 16. Since 
the correlation functions are by definition (weighted) aver- 
ages over the full multipole range, the reduced sensitivity 
to this type of non-Gaussian structure is certainly not un- 
expected. 
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